The cosine law at the atomic scale: Toward realistic simulations of Knudsen diffusion 
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We propose to revisit the diffusion of atoms in the Knudsen regime in terms of a complex dynam- 
ical reflection process. By means of molecular dynamics simulation we emphasize the asymptotic 
nature of the cosine law of reflection at the atomic scale, and carefully analyze the resulting strong 
correlations in the reflection events. A dynamical interpretation of the accomodation coefficient 
associated to the slip at the wall interface is also proposed. Finally, we show that the first two 
moments of the stochastic process of reflection non uniformly depend on the incident angle. 
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A century ago, Knudsen proposed a model describing 
the diffusion of dilute gases through cylindrical pores Q . 
The mean free path of the molecules being larger than the 
pore size, the transport properties are essentially driven 
by the collisions with the wall. In the Knudsen approach, 
based on the kinetic theory of gases, all rebounds were 
assumed to be governed by the diffusive Lambert's co- 
sine law of reflection [2]. Later, the model was improved 
by Smoluchowski who proposed, using Maxwell's theory, 
that only a fraction / of rebounds were diffuse, while the 
remaining 1 — / were specular [3|. The physical situ- 
ations addressed by these authors gain renewed interest 
due to the emerging applications in gas separation [4] and 
catalysis by means of new porous materials [5]. In order 
to achieve efficient processes a thorough understanding 
of the mechanisms involved in transport of gas through 
porous membrane is required. This problem has essen- 
tially been addressed by two distinct types of numerical 
simulations: billiard-like simulations (BLS) and molecu- 
lar dynamics (MD). In the first, the gas- wall interaction 
is introduced through ad hoc laws of reflection, disregard- 
ing microscopic mechanisms. Therefore, BLS is more ap- 
propriate to describe the impact of the pore geometry on 
the transport properties at a macroscopic scale [gI, (tI, (sj . 
In MD microscopic interactions are properly considered, 
the cost being a limitation in the spatial and time scales 
investigated by the simulations [1, [lH • 

To combine the advantages of both kinds of simula- 
tions, one could extract from MD the realistic reflection 
law that would be finally incorporated in BLS. A first 
step in this direction has recently been made by Arya et 
al. Using MD, these authors quantified the rela- 

tionship between the phenomenological coefficient / and 
the parameters of the Lennard- Jones potential describ- 
ing gas- wall interactions, namely the wall structure and 
the interaction energy. Nonetheless, a realistic way of 
introducing / in BLS has not been achieved yet. 

In this letter we propose to fill the gap between MD and 
BLS simulations. First, we show the robustness of the 
cosine law by successfully testing it against MD results. 
Then, the velocity distribution functions are analytically 
derived and successfully compared to numerical results. 



We re-examine the cosine law of reflections emphasizing 
its asymptotic nature. Through a careful analysis of the 
correlations between consecutive reflection events, we es- 
tablish a relation between / and the range of the dynam- 
ical correlations. Finally, we present an analysis of the 
reflection law in terms of a complex stochastic process. 




FIG. 1: Snapshot of the simulated system where a rebound 
event is schematically depicted. The dotted line shows the 
border between the ballistic and interacting regions. 



We use MD simulations to describe the trajectory of a 
particle through a two-dimensional slit pore, taking into 
account the atomistic structure of the walls (see Fig. [T]). 
The equations of motion are integrated through a stan- 
dard Verlet Algorithm [13] with periodic boundary con- 
ditions in both directions. The two different interaction 
potentials, between atoms composing the wall and be- 
tween the diffusing particle and the atoms, are of the 
Lennard- Jones type: V{r) = 4e{{(7/ry^ — (<j/r)^]. We 
use the shifted force potential to guard against energy 
and force discontinuities at the cut-off radius Tc = 2.5cr 
[Tsj . For both potentials a = 1, for interactions between 
the wall atoms e = 1 while e takes values ranging from 
0.005 to 0.5 for the particle-atoms potential. The mass of 
the atoms is fixed to m = 1. As a consequence, lengths, 
energies and masses are respectively expressed in units 
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The walls are composed of 141 atoms arranged in a 
triangular structure (see Fig. [1]), the height of the pore 
is fixed to 6<j. Since the potential is of finite range, we can 
easily define the ballistic region, in the central part of the 
pore, where the diffusing particle does not interact with 
the walls. To characterize the "rebounds" we therefore 
record the particle positions and velocities at the borders 
of this region (see Fig. [T]). The Boltzmann constant is 
put to unity and the temperature, expressed in units of 
energy, is fixed to T = 0.15, a value for which the solid 
phase of the wall is stable. Once the system is in thermal 
equilibrium we perform simulations in the microcanonical 
ensemble. A typical run consists in 3 x 10^ time steps 
{At = 10~^) during which we roughly compute from 10^ 
to 4 X 10^ rebounds on the walls depending on the value 
of e. 




FIG. 2: Probability density function of the rebound angle. 
Simulation data (circles) obtained for e = 0.1 are in complete 
agreement with the expected cosine law {[T]) (full line). The 
same agreement is observed whatever the value of e. 

The probability density function of the refiection an- 
gles p{0) obtained numerically through an averaging over 
2 X 10^ bounces and for e = 0.1 is shown in Fig. [2l Our 
numerical results confirm the expected cosine law of re- 
fiection: 

p{e)de = i cos OdO, p(sin 0)d sin = ^d sin (1) 

The same remarkable agreement is observed for all the 
values of e used in our simulations, showing the robust- 
ness of the cosine law. It is worth noting that relation ([T]) 
represents an asymptotic law (time t ^ oo) which does 
not preclude correlations between successive angles of re- 
fiection. One can even obtain such "ergodic" behavior 
with a pure specular law of refiection, for a non-trivial ge- 
ometry of the surrounding walls. Indeed, in the limiting 
case of billiard problems, which constitute paradigmatic 
examples of Hamiltonian chaos, one generically obtains 
the cosine law of refiection whereas incident and refiected 
angles are specularly related [14]. The role of short-time 



correlations between angles of refiection, which cannot 
be neglected, will be addressed in detail in the last part 
of the paper. 

We use the cosine law as the starting point to derive 
analytically the other distribution functions describing 
the motion of the diffusing particle in the ballistic zone. 
As expected, our simulations show that the velocity com- 
ponent parallel to the surface, Vx = vs'mO, is a Gaussian 
random variable with variance given by the fixed tem- 
perature (for simplicity, the mass of the particle is 1): 
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The velocity can be described by two sets of random vari- 
ables: (vx^ Vy) or ('u, sin^) whose distributions are related 

by 
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We numerically checked that Vx and Vy (or, equiva- 
lently, v and sin^) are independent random variables: 
p{vx,Vy) = px{vx)py{vy) (p(v, sin i9) = py{v)pe{s\nO)). 
Thus, the distribution py can be obtained from 
the probability density of the modulus Pv('^) = 
J^-j^ d(sin^) p('U, sin ^). From the normalization con- 
straint dv py = 1 one, indeed, obtains: 



d(sin6>) / dv p('u cos 6>) exp 

cos 



sin^ i 



2T 



= (4) 

where we used the Maxwell-Boltzmann distribution 
([2j). The solution of (|4]) is given by p{vcosO) = 
{v cos6/T) exp [—v'^ cos^ 0/2T] . Finally, we can assume 
for Vy the following distribution 



PyM = ^exp 



2T 



(5) 



Due to the non-symmetric interactions experienced by 
the particle in the ^/-direction, the corresponding velocity 
distribution ([5]) differs qualitatively from ([2]). Recall that 
([5]) is fully compatible with the diffuse law of refiection 
([1]). A similar expression of ([5]) had been phenomenolog- 
ically obtained by Arya et al. [12,]. 

Knowing px and p^, we are now able to write the ex- 
plicit form of py'. 
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As expected from (|5|), it is clear that the modulus v 
does not obey the standard two-dimensional Maxwell 
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distribution \y?\. Expression (|6]) is successfully tested 
against the simulation results as shown in Fig. [3l Here 
again the agreement is quite remarkable whatever the 
values of e. For comparison the 2D Maxwell distribu- 
tion is also shown in Fig. [3l Note that the mean value 
{v) = y^ST/tt deduced from (|6|) is around 20 % smaller 
than the Maxwell mean value (vm) = We want 

to draw attention to the fact that distribution (|6|) should 
be used in BLS of Knudsen diffusion. Indeed, the com- 
monly used Maxwell distribution leads to an underesti- 
mate of the diffusion constant. 



1 ^ r~ 

° Simulation 

— P(v) 

— Maxwell 




0^ 1 
Velocity 



1^ 



FIG. 3: Probability density function of the velocity modulus. 
Circles are simulation data, the full and dotted curves corre- 
spond respectively to pv{v) given in (|6)) and to the Maxwell 
distribution. 

We now turn to a thorough analysis of correlations be- 
tween successive reflected angles in order to understand 
how the cosine law of reflection is asymptoticaly attained. 
We thus define the following function: 



C{n) = 1 



(7rV2)-4 



(7) 



C{n) gives a measure of the correlation between two an- 
gles of reflection separated by n rebounds of the particle 
on the walls. By construction, C{n) = 1 for a pure spec- 
ular regime of reflections and C{n) = in the diffusive 
regime. The evolution of C{n) for three values of e is 
shown on Fig. [H We observe that the expected conver- 
gence toward the diffusive limit as n increases is observed. 
We numerically verified that the mean resident time of 
the particle inside the interacting region increased with 
e. As a consequence, one can see on Fig. 2] that the rate 
of convergence to C = increases with larger e. But it is 
worth noting that C{n) exhibits non-zero values for a sig- 
nificant number of bounces. The correlation functions ex- 
hibit an exponential behavior C{n) = exp{—n/nc) where 
the characteristic number of rebounds Uc gives the range 
of the correlation: two angles of reflection separated by 
a few Tie rebounds are uncorrelated. For the correla- 
tion functions shown in Fig. HJ Uc varies from 0.23, for 
the largest value of e, to 2.92, for the smallest. Such a 
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FIG. 4: Correlation function C{n) represented for different 
values of e. Circles, diamonds and squares correspond to e = 
0.005, 0.02 and 0.1. Inset: C{1) is plotted as a function of e 
respectively. 



large number of "non-diffusive" (or "quasi-specular" ) re- 
bounds may drastically alter the transport process. This 
is what had been phenomenologically taken into account 
in the Smoluchowsky model [3] by introducing fractions / 
and 1 — / of atoms having diffuse and specular rebounds 
respectively. This fraction / is called the tangential mo- 
mentum accommodation coefficient and can be related 
to the slip coefficient at the wall [1, [lij . Our approach 
gives an insight into the microscopic dynamics which gen- 
erates the observed macroscopic behavior and provides a 
new interpretation of the phenomenological coefficient. 
In order to mimic the observed correlations, one can in- 
deed imagine a more realistic BLS where Uc rebounds 
are specular, after which the next rebound is diffuse and 
followed by another Uc specular rebounds. Thus, the ac- 
commodation coefficient can be viewed as the inverse of 
the characteristic number of rebounds: / ~ l/nc- 

Let us now focus on C{1) which quantifies the loss of 
memory after a single bounce. The inset in Fig. 2] shows 
the evolution of this "instantaneous" correlation function 
with the strength of the potential. Even for the small- 
est value of e the reflection is never truly specular (in 
this case C{1) = 1) due to the atomic structure of the 
wall. As expected, for increasing values of e C{1) van- 
ishes, indicating a complete randomization of the reflec- 
tion events. This is clearly illustrated by looking at the 
motion of the particle for e = 0.1 where the absence of 
correlations between incident and reflected angles is evi- 
dent (see EPAPS 0, video 1). Note that the diffusive 
behavior stems from the fact that for almost all rebounds 
the particle experiences, in the interacting zone, multi- 
ple collisions with wall atoms. This is no longer the case 
for 6 = 0.01 (see EPAPS flG*], video 2) where a unique 
collision event generates the rebound. 

Since the reflected angle is never completly diffusive or 
specular but something in between, the important quan- 
tity to know is the conditional probability p{0r\0i) to re- 



4 




-jt/2 -jt/4 „ k/4 nil 



FIG. 5: First and second (inset) moments of the distributions 
pe^{AO). The full lines correspond to the purely diffusive 
regime. Circles, diamonds and triangles correspond respec- 
tively to e = 0.005, 0.05 and 0.5. 



proximations in the studies of Knudsen diflFusion: which 
(z) assume Maxwell velocity distribution functions, (ii) 
neglect the correlations between consecutive rebounds. 
We gave the correct distributions, Eq. ([5]) and (|6|), ac- 
tually involved in the Knudsen regime. By defining an 
appropriate correlation function ([7j), we carefully ana- 
lyzed the strong correlations occurring in the dynamical 
process of reflection. We revealed a characteristic cor- 
relation rebound number ric that has been linked to the 
commonly used accommodation coefficient /. Finally, we 
performed the complete characterization of the stochastic 
rebound process and emphasized the non-trivial depen- 
dence of the conditional probability p{0r\0i) on the inci- 
dent angle. This work should pave the way for new stud- 
ies of Knudsen transport using improved BLS schemes as 
proposed in this letter. 

We gratefully acknowledge fruitful discussions with A. 
ten Bosch. We also wish to thank G. Batrouni, A. ten 
Bosch and Sebastien Tanzilli for their careful reading of 
the manuscript. 



bound with an angle Or, the incident angle Oi being fixed. 
If we want to focus on the deviation from the specular 
regime, it is convenient to write Or = —Oi-\-AO and there- 
fore express the conditional probability as a set of distri- 
butions pe,{AO). We computed {AO{Oi)) and {AO{Oi)^), 
the first and second moments of these distributions, and 
plot it for different values of e in FigO In the case of a 
purely diffusive reflection we expect: 

(A^(^,)) = 0, (8) 
{AO{Oif) = ^2 + ^2/4_2 (9) 

As shown in Fig. [5l the diffusive behavior is recovered 
for the most attractive potential (e = 0.5) considered in 
this study. The opposite case of purely specular reflec- 
tions, characterized by {AO{Oi)) = and {AO^Oif) = 0, 
is never reached. Indeed, as discussed above when in- 
terpreting the fact that C(l) 7^ 1 even for e ^ 0, the 
wall structure precludes the emergence of a pure specu- 
lar regime. In the intermediate regimes (e = 0.05, 0.005), 
one can see that both moments non uniformly depend 
on the incident angle, the non-generic behavior being 
more important for the largest incident angles. Thus, we 
demonstrate here that the correlations between two suc- 
cessive rebounds are quantitatively sensitive to the value 
take by the first rebound. We are currently developing an 
analytical description of such a complex behavior in the 
choice of the reflected angle in an improved BLS scheme. 

In this letter, we have proposed new insight into the 
so-called cosine law of reflection by means of MD simula- 
tions. We confirmed its validity at the atomic scale what- 
ever the strength of the gas-wall interaction, but only as 
an asymptotic law. We examined two crude "ergodic" ap- 
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